/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */


/*! \file D_c_.h
    \brief Multidimensional generalization of centered finite difference operators on hypercube.

	Multidimensional generalization of the classical centered difference operator
	\f[
			\delta_{\hat{x}} u_{i}=\frac{u_{i+1}-u_{i-1}}{2h}
	\f]
	The geometry is the unit cube \f$[0,1]^N\f$, where N is given by data::Dim
*/

#if !defined(_CenteredDifference_)
#define _CenteredDifference_
#include <src/math/BoundaryCondition/classicalgraphdefinition.h>
#include <src/math/Operators/FiniteDifferenceOperators.h>
#include <src/algorithm/Evolution_Algorithm/CFL_traits.h>
#include <src/algorithm/Evolution_Algorithm/Predictor_traits.h>

/*! 
    \brief One dimensional specialization of centered difference operators.

	\f[
			\delta_{\hat{x}} u_i=\frac{u_{i+1}-u_{i-1}}{h}
	\f]
 */
template <class data, class V = smart_ptr<data>>
class D_c_  
	: public FiniteDifferenceOperator<data,data,V>
{
public :
	OS_STATIC_CHECK(data::Dim==1);
		enum:std::size_t												{Which_Dim = 1};
		enum:std::size_t												{Order = 1};
		typedef OTSCentered												Centered_Type;
		typedef typename data::Boundary_type							Boundary_type;
		
		typedef FiniteDifferenceOperator								Base;
		typedef  typename Base::Base_Cast								Base_Cast;
		typedef D_c_													type;

//! return " D_c u".
		static std::string get_name(){return " D_c u ";};
		D_c_( ) : operator_interface_(*this,make_2<base_value_type>(0.,1.)) {};
		D_c_(const D_c_ & data) : operator_interface_(*this,make_2<base_value_type>(0.,1.)),  Base(data) {};
		D_c_(const weight_set_type & data) :  operator_interface_(*this,make_2<base_value_type>(0.,1.)), Base(data) {};

		D_c_(const Base_Cast &first) : operator_interface_(*this), Base(first)
		{};
		virtual ~D_c_(){};
		D_c_ & operator = (const D_c_ & from) {
			Base::operator =(from);
			return *this;
		};
		virtual void update() {
			operator_interface_.set_size(get_state_variable()->size(Which_Dim));
		};
		//! Optimized code for centered finite difference operator.
		virtual result_type operator()(typename argument_type from) const //optimization code
		{
			return ((from<<1) - (from>>1) )*(.5*from.size());
		};
		//! Factory method for Algebra.
		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper_arg<type,weight_set_type>(&weight_set_);
		};
private :
		operator_interface<D_c_,Order,Centered_Type,Boundary_type> operator_interface_;
};
/*! \brief Bind operator of centered difference operator D_c_

	This is the time Bind operator of backward difference operator D_c__traits, that is D_c__t(t,u)=D_c_(u)
*/
template <class data, typename Time = OS_double >
struct D_c__t : public BindTimeOperator<D_c_<data> ,Time > {

	typedef D_c__t															time_operator_type;
	typedef BindTimeOperator<time_operator_type ,Time >						Base;
	D_c__t(){};
};

/*! 
    \brief multidimensional finite difference centered operators.

	This is the centered discrete gradient vector. More precisely : 
	\f[
			\nabla_{\hat{x}} u_{i_1,...,i_D}=\{ \frac{u_{...,i_j+1,..}-u_{...,i_j-1,..}}{2h_j} \}_{j=1,..,D}
	\f]
*/
template <class A, class R, std::size_t D = A::Dim, class V = smart_ptr<A> >
class Grad_c_ : public OS_Generic_Algebra<A,R,V>
	, public virtual graph_base<A>
 {
 public :
		enum:std::size_t												{Dim = D};
		OS_STATIC_CHECK(Dim==2); // lazyness...
		enum:std::size_t												{Order = 1};
		typedef OTSCentered												Centered_Type;
		typedef typename argument_base_type::value_type					argument_value_type;
		typedef typename argument_value_type::Boundary_type				Boundary_type;
		
		typedef OS_Generic_Algebra										Base;
		typedef  typename Base::Base_Cast								Base_Cast;
		typedef Grad_c_													type;

//! return " D_c u".
		static std::string get_name(){return " D_c u ";};
		type( ) {};
		type(const type & from) : Base(from) {};
		type(const Base_Cast &first) : Base(first)
		{};
		virtual ~Grad_c_(){};
		type & operator = (const type & from) {
			Base::operator =(from);
			return *this;
		};
		virtual void update() {
		};
		virtual result_type operator()(argument_type from) const 
		{
			typedef argument_base_type::const_iterator const_iterator;
			typedef argument_base_type::Base Base_argument_type;
			typedef result_type::value_type::iterator iterator;
			typedef result_type::indice_list_type indice_list_type;
			if (from.capacity() == 0) return result_type();
			indice_list_type ind(Dim+1);
			std::copy(from.size_list().begin(), from.size_list().end(),ind.begin()); 
			ind[Dim] = Dim;
			result_type result(ind);
			argument_base_type::size_type size1 = result.size(1);
			argument_base_type::size_type size2 = result.size(2);
				OS_Array_Copy(result[0], ((from<<1) - (from>>1) )*(.5*size1));
				const_iterator beg = from.begin();
				const_iterator end = from.end();
				iterator beg_r = result[1].begin();
				for (;beg != end; ++beg, ++beg_r) {
					*beg_r = ((*beg<<1) - (*beg>>1) )*(.5*size2);
				};
			return result;
		};
		//! Factory method for Algebra.
		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper_arg<type,type>(this);
		};
};
/*! 
    \brief multidimensional finite difference centered operators.

	This is the centered discrete divergence vector. More precisely : 
	\f[
			\nabla_{\hat{x}} \cdot u_{i_1,...,i_D}=\sum_{j=1,..,D} \frac{u_{...,i_j+1,..}-u_{...,i_j-1,..}}{2h_j} 
	\f]
*/
template <class A, class R = typename A::value_type, std::size_t D = R::Dim, class V = smart_ptr<R> >
class Div_c_ : public OS_Generic_Algebra<A, R,V>
	, public virtual graph_base<A>

 {
 public :
		enum:std::size_t												{Dim = D};
		OS_STATIC_CHECK(Dim==2); // lazyness...
		enum:std::size_t												{Order = 1};
		typedef OTSCentered												Centered_Type;
		typedef typename A::Boundary_type								Boundary_type;
		
		typedef OS_Generic_Algebra										Base;
		typedef  typename Base::Base_Cast								Base_Cast;
		typedef Div_c_													type;

//! return " Grad_c_ . u".
		static std::string get_name(){return " Grad_c_ . u ";};
		type( ) {};
		type(const type & from) : Base(from) {};
		type(const Base_Cast &first) : Base(first)
		{};
		virtual ~Div_c_(){};
		type & operator = (const type & from) {
			Base::operator =(from);
			return *this;
		};
		virtual void update() {};
		//! Optimized code for centered finite difference operator.
		virtual result_type operator()(typename argument_type from) const 
		{
			typedef argument_base_type::const_iterator const_iterator;
			typedef result_type::iterator iterator;
			if (from.capacity() == 0) return result_type();

			argument_base_type::size_type size1 = from.size(1);
			argument_base_type::size_type size2 = from.size(2);

			result_type result(size1,size2);
			const_iterator beg = from.begin(), end = from.end();
			OS_Array_Copy(result, ((*beg<<1) - (*beg>>1) )*(.5*size1));	//partial_1 f^1


			for (++beg;beg != end; ++beg) {
				typedef result_type::iterator iterator;
				iterator beg_1 = result.begin(), end_1 = result.end();
				typedef argument_base_type::value_type::const_iterator const_iterator;
				const_iterator beg_beg = (*beg).begin();
				for (;beg_1 != end_1; ++beg_1, ++beg_beg) {
					(*beg_1) += ((*beg_beg<<1) - (*beg_beg>>1) )*(.5*size2);   //partial_2 f^2 and so on...
				};
			};

			return result;
		};
		//! Factory method for Algebra.
		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper_arg<type,type>(this);
		};
};

#endif